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Abstract 

A model is proposed to study the process of hypoxia-induced angiogenesis in 
cancer cells. The model accounts for the role played by the vascular endothelial 
growth factor (VEGF)-A in regulating the oxygen intake. VEGF-A is dynami- 
cally controlled by the HIF-la concentration. If not degraded, HIF-la can bind 
to the subunit termed HIF-1/3 and so experience translocation to the nucleus, 
to exert its proper transcriptional activity. The delicate balance between these 
opposing tendencies translates into the emergence of distinct macroscopic be- 
haviors in terms of the associated molecular concentrations that we here trace 
back to normoxia, hypoxia and death regimes. These aspects are firstly ana- 
lyzed with reference to the ideal mean-field scenario. Stochastic fiuctuations 
are also briefly discussed and shown to seed a cooperative interaction among 
cellular units, competing for the same oxygen reservoir. 
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1. Introduction 



Hypoxia (O2 tension below 2.5 mmHg) is a hallmark of several types of solid 
tumor [ll| . Growing clinical evidence postulates a correlation between hypoxia 
and cancer aggressiveness and metastasis 0] ■ Hypoxia may underlie resistance 
to radiotherapy and can modify the cancer phenotype through gene regulation. 
A major cellular response to hypoxia is the stabilization of hypoxia inducible 
factor 1 (HIF-1), a transcription factor that controls cancer progression, and 
represents a potential target for therapy HIF-1 is a heterodimer composed 
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of a oxygen-dependent HIF-la and a constitutively expressed subunit termed 
HIF-1/3. The transcriptional response to HIF-1 a varies from general effects, 
such as the up regulation of anaerobic respiration (virtually in all tumors), 
to the more tissue-specific effects such as angiogenesis [llj. The tumor-driven 
angiogenesis is therefore a process pivotal in tumor growth and progression, 
and is mainly related to the HIF-1 dependent secretion of the angiogenic factor 
Vascular Endotclial Growth Factor -A (VEGF-A) [l^]- Low O2 tension may 
shape the cancer phenotype. According to a recently proposed model for cancer 
progression the tumor milieu exerts Darwinian selection in favor of cancer 
cells. The basis for this selection can be two-fold: firstly, the milieu may activate 
tumor-specific pathways. A second basis for selection may reside in the ability 
of cells to protect the intracellular compartment from a hypoxic milieu. Given 
the important role of hypoxia in cancer progression targeting the cellular 
responses to hypoxia, including the triggering of VEGF-A secretion and the 
ensuing angiogenesis, may form an alternative to, or a component of, the current 
chemothcrapcutical practice. 

In this paper, a model is proposed to resolve the complex dynamical interplay 
between the molecular chemical species that participate to the aforementioned 
process. In particular, we will start by exploring the mean-field scenario stem- 
ming from the hypothesized networks of interaction. This formally amounts to 
operate in the limit for infinite system size, the species concentrations obeying 
to a closed set of ordinary differential equations. We target our analysis to the 
asymptotic (equilibrium) regime highlighting the emergence of different attrac- 
tive (stable) states of the single cell dynamics, which appear to be selectively 
chosen depending on the chemical parameters involved. We then turn to con- 
sider the evolution of three neighbors cells which are effectively coupled via the 
external environment, and adjust consequently their interior dynamics. The role 
of fluctuations stemming from the discreteness of the system under scrutiny, is 
shortly addressed resorting to direct stochastic simulations. Depending on the 
selected parameters, normoxic condition can develop as an emergent dynami- 
cal equilibrium. Allowing for punctual mutations, so to mimic tumor derive, 
can alter the stability of the system, and eventually result in a competition be- 
tween distinct cell populations. Indeed, as we shall demonstrate, cooperative 
mechanisms arise being mediated by the stochastic component of the dynamics. 
Interestingly, adjacent cells can also profit from the gained ability of an individ- 
ual cell unit to survive under hypoxic condition, an effect which fades off in the 
continuum limit. 

The paper is organized as follows. Next section is devoted to introducing 
the key ingredients of the model, here formulated as an ensemble of chemical 
equations. Then the mean- field version of the model is studied, which in turn 
corresponds to work in the infinite system size limit. The focus is on a single 
cell: different dynamical regimes are identified corresponding to distinct choices 
of the chemical parameters. Stochastic simulations are also performed to test 
the adequacy of the theory predictions and challenge in silico the role played by 
finite size corrections. Then in Sec. [5] we turn to discussing a generalization of 
the proposed model where three cells are made to interact and compete for the 
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oxygen amount. Finally, in Sec. ^ we sum up and draw our conclusions. 

2. A chemical model of hypoxia 

Consider now an individual cell and assume the membrane barrier to divide 
between interior and exterior regions. Selected molecules can migrate through 
the semi-permeable membrane, being hence transported from the outside to 
the inside and viceversa. Aiming at elaborating a sound description of the 
hypoxia cycle [l^, we select a set of candidate molecules which define a close 
ensemble of interacting elements, according to our schematic representation of 
the phenomenon under scrutiny. Oxygen molecules (O2) populate the external 
milieu. They diffuse and eventually reach the membrane where they happen 
to combine with the hydroxylases (W), which are consequently turned into 
an active phase, hereafter labeled Wq0. This process is encapsulated into the 
following chemical equation: 

O2 + W ^ Wa + Eo (1) 

where a is the reaction rate and Eq stands for the empty case which the oxygen 
leaves behind when it gets absorbed by the W. The label refers to the chemicals 
which are populating the outside region, while the molecules confined inside the 
cell wall are targeted with the subscript i. By introducing the concept of empty 
spaces (Eq outside, and E, inside) we impose that the total number of elements 
(including the empties) is a constant of the dynamics Q ■ Active hydroxylases 
Wa interfere with the HIF-a, yielding to the degradation via the proteasome 
pathway: 

HIF-a + W„ -^W + E, (2) 

Once the HIF-a is being targeted to degradation, the hydroxylases go back 
to the primitive configuration, waiting for the next oxygen to drive the transition 
into its active state. The degraded HIF-a molecule is replaced by the empty 
inner element E^. 

Moreover, the HIF-a can occasionally meet the HIF-/3 subunit. This en- 
counter necessarily takes place within the cell and gives rise to the VEGFj 
macromoleculcH. The corresponding chemical equation reads: 

HIF-a 4- HIF-/3 VEGF, + E, (3) 

VEGFj can eventually abandon the cell interior to occupy an external va- 
cancy Eq, which is hence mutated into a VEGFq element. In formulae: 



^Hydroxylases arc intracellular enzymes, organized in three families respectively PHDl, 
PHD2, PHD3. Hydroxylases hence populate the internal flow and eventually combine with 
the oxygen that have penetrated the cell. As formalized in the following, we shall here focus 
on a more simple picture, and imagine that the encounters between hydroxylases and oxygen 
occur at cell boundary. In such a way, hydroxylases can be ideally pictured as antennae, 
localized on the membrane wall, chasing for the (outer) oxygen molecules. 

^From hereon we will simply refer to VEGF-A as to VEGF. 
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VEGF, + Eo ^ VEGFo + E, 



(4) 



Furthermore, the VEGFq participates to the network of extracellular reac- 
tions which underlay the angiogenesis process. In practice, following a complex 
cascade of nested chemical reactions, the VEGFq attracts the oxygen O2, which 
is hence imagined to replace an (external) empty case. This materializes in turn 
into a simplified vision of a complex dynamic pathway, which has the sole scope 
of providing a self-consistent entry to the inspected process. The corresponding 
chemical equation takes the form: 

VEGFo + Eo ^ VEGFq + O2 (5) 

Notice that the VEGFq molecule is also a product of the reaction, a working 
hypothesis that we motivate with the empirical observation that VEGFq can 
attract several oxygen units. This process will eventually come to a stop, as 
follows VEGF degradation which occurs both inside and outside the cells: 

VEGFq Eq (6) 

VEGF, E, (7) 

Finally, we should also account for the constitutive generation of HIF-a and 
HIF-/3: 

E, ^ HIF-a (8) 
E, HIF-/3 (9) 

The scenario outlined above is depicted in Fig. [U where the main molecular 
actors are also specifiecH. It should be noted that the above model conserves 
the number of molecules which totals in A^, a fact that is ultimately related to 
the presence of the empties species E. 

We should emphasize that the proposed model is intimately stochastic, in its 
original chemical formulation. The inherent stochasticity emerges as an effect of 
the discreteness of the medium: finite size corrections, also termed demographic 
noise, result in a endogenous perturbation which can significantly impact the 
dynamics of the system as compared to the idealized continuum picture, which 
formally applies in the limit of infinite constituents. To respect the stochastic 
nature of the problem and simulate it in silico, one can resort to the celebrated 
Gillespie algorithm , a method which makes it possible to exactly monitor the 
time evolution of the discrete population of mutually interacting entities. The 



■^Importantly, the oxygen is initially imagined to fill the external reservoir, and solely re- 
integrated via the VEGF pathway. This is clearly a simplifying assumption, as in general one 
should also account for the oxygen supply via blood fiux. On the other hand, we are here 
interested in the ability of the cell to self-organize when exposed to stressing condition, as it 
is the lack of oxygen income, a working hypothesis which justifies the assumed scenario. 
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Figure 1: Cartoon for a single coll dynamics as follows the proposed scheme. Oxygen molecules 
(circles laying outside the cell) are captured by the hydroxylases sitting on the membrane, 
and so trigger the HIF-a (dark circles, red online, stored inside the membrane) elimination 
pathways (dashed circle). At the same time, HIF-a can combine with HIF-/3 (light circles, 
orange online, stored inside the membrane), yielding to the (barrel like, in the picture) VEGF^. 
This latter can exit the cell wall, so ensuing in VEGFq and drive the incoming oxygen flux. 



scheme captures in fact the probabihstic nature of the microscopic couphngs 
and enables one to resolve the contribution associated to finite size fluctuations. 
As we shall demonstrate, stochastic fluctuations can eventually materialize in 
complex dynamics, which call for sound biological interpretation. 

As anticipated, and opposed to the probabilistic approach, one can invoke 
the deterministic approximation, which in turn corresponds to assuming an in- 
finite population amount and consequently disregarding finite size correlations. 
This defines the mean-field, continuum level of description, to which next sec- 
tions are entirely devoted. More specifically, we will start by considering the 
case of a single cell, which we will discuss with reference to its mean field ana- 
logue. Notice that the deterministic formulation results from a straightforward 
derivation, that moves from the underlying master equation for the stochastic 
process, as defined by the above set of chemical equations. In the following, we 
will allude to such procedure, recalling its general philosophy, and just provide 
the end result for the case under inspection, without insisting on the techni- 
calities. Working within the continuum representation we will then elaborate 
on the asymptotic dynamics of the recovered system of ordinary differential 
equations, by identifying distinct fixed points and discussing their associated 
stability characteristics. The inherent stochastic effects stemming from finite N 
corrections, will be relegated to section 5. 
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3. Master equation 

At time t, the status of the chemical system, as ruled by the above reactions, 
is completely specified once the species amount are assigned. These are discrete 
entries, one for each species belonging to the ensemble of interacting population, 
which define the scalar components of the status vector n. Recall that the 
system is intrinsically stochastic, its evolution being a suite of events which 
realizes with a prescribed probability. In mathematical terms, this fact implies 
dealing with the probability that the system is seen in the state n at time t, 
hereafter denoted P(n, t). This latter quantity adjusts in time and obeys to the 
master equation: 

|p(n, t) = Yl [nn\n')P{n', t) ~ T(n'|n)P(n, t)\ (10) 

n' 

where T(n'|n) denotes the transition rates between two adjacent states, from 
n to n', compatible with the chemical constraints imposed by equations ([1])- 
([9]). For the sake of compactness, we omit here to list the specific forms of the 
transition rates relative to the model under analysis. 

The partial differential equation (jlOp is an exact formulation, which can be 
in principle analyzed to track the system dynamics. In practice however, it is 
hard to analytically handle and further progress has to rely on direct numerical 
investigations. 

The Gillespie algorithm represents in particular a viable strategy to simulate 
the time evolution of the probability distribution P(n, t). The method sets up a 
Monte Carlo procedure to determine the next reaction that is selected to occur, 
among all possible ones, and the time r, when the event takes place. Both 
selections reflect the assigned transition probabilities, which scales proportional 
to the number of substrate molecules times the associated reaction constants. 
Then, based on the selected reaction, the update of the molecule count follows, 
and time t is replaced by its value increased by the stochastic time interval t, 
namely i — > t + r. 

Alternatively, the master equation (|10p can be studied via perturbative 
methods, aimed at extracting both the idealized continuum picture and the 
successive finite size corrections. The van Kampen large TV expansion is in this 
respect a widely adopted technique, extensively characterized in the literature 
■ In the following section, we shall operate within the van Kampen frame- 
work to derive a closed system of ordinary differential equations for the coupled 
evolution of the concentrations of chemical system ([I])-® in the limit of large 
(formally infinite) N . 

4. Mean-field solution of the single cell model 

Let us focus on a single cell, as pictorially represented in Fig. [1] Consider for 
instance the O2 species. According to van Kampen [l^ . one can put forward 
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the following ansatz for the associated mean-field (continuous) concentration 
here labeled O: 



where the extensive quantity is discrete: it is one of the components of 
the vector n and stands for the number of oxygen molecules. ^ is instead 
a stochastic variable and contributes with an additive terms which scales as 
so materializing in a perturbation to the idealized mean-field dynamics. 
Stochastic effects can be hence safely neglected when performing the continuum 
limit 0, [ill , the relic concentration obeying to a deterministic set of coupled 
differential equations. We recall however that beyond the idealized scenario, 
finite size corrections do matter when real systems are concerned, the intimate 
discreteness being a crucial ingredient of the dynamics. We shall return on this 
issue in Section 5, to highlight the important impact of such an endogenous 
source of stochastic perturbation with reference to the scrutinized model. 

For simplicity, by extending the above reasoning to the other involved species, 
wc will here denote by W, Wa, Vo, Vi, £o, £i, Ha and Hp the average (mean- 
field) concentration of W, Wa, VEGFq, VEGFj, Eq, Ei, HIF-a and HIF-;3 re- 
spectively. According to this notation, the mean-field equations take the forrrQ: 

j^O = -aOW + fVo£o (11) 

= ~aOW + bWaHa (12) 

jWa = aOW-bWaHa (13) 

jHa = -bWaHa - cHaH/i + ga£^ (14) 

^■H/j = -cHaHp + g/)£i (15) 

jV, ^ cHaHfj - eV,£o - dV, (16) 

^Vo = eV,£o-dVo (17) 
dt 

j^£, = bWaHa + cHaHi3 + e£oV, + dV^- {ga+gp)£, (18) 
with the conservation law 

O + W + Vo + V^+£o + Wa+£^+Ha+Hp ^1 (19) 



*In deriving eqs. II18I I correlations have been neglected, a legitimate choice in the mean-field 
approximation and whose validity has been tested versus numerical simulation in a number 
of models that follow the same conceptual scheme 0, . 
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Notice that the same result could be recovered without the use of the van 
Kampen's expansion. Indeed, multiplying both sides of the master equation by 
a scalar component of the status vector n, summing over all n, shifting some of 
the involved sums by ±1 and neglecting the correlations (which is legitimate in 
the N infinite limit) one recovers the same mean field system. 

From Eq. and ((131), it straightforwardly follows that + W stays 

constant: This invariant quantity is hereafter set equal to no and specified by 
the initial condition. To investigate the dynamics of the above system (fTTj) - 
P^ . and in particular elaborate on its asymptotic evolution, one can look for 
stationary solutions. These are the fixed points of the dynamics, which can 
be explicitly found by setting to zero the derivatives in ([TT|) - and solving 
the obtained algebraic system, where the (final) concentration amounts enter as 
unknowns. 

4.I. Stationary states 

Solving for the fixed points returns three possible solutions. These latter are 
here characterized and labeled via the progressive indexes i = 1,2,3, associated 
to the predicted asymptotic concentrations. Two solutions can be readily ob- 
tained via direct inspection of the above equations, and correspond to opposite 
tendencies that we will hereafter call normoxia and death respectively. 

Normoxia conditions are assumed to correspond to the asymptotic solution: 



HIF-a is turned off and so arc the inner/outer concentrations of VEGF. The 
HIF-a degradation mechanism hence prevails over the alternative pathway that 
yields VEGF production, and which could in principle contribute to enhance 
the oxygen supply. The amount of O2, here quantified by Oi, is on the contrary 
sufficient to guarantee the correct cell functioning. Hydroxylases arc frozen in 
their active state Wa and consequently stimulate the repression of the HIF-a 
quotgd. From the conservation relations we immediately deduce the following 
relation: 



The outer cases are shared between oxygen and empties, as exemplified by the 
latter equation, in a proportion that cannot be predicted a priori. 



^In the transient dynamics, HIF— o molecules are indeed present and translate into VEGF 
entities, upon combination with HIF— /3. VEGF molecules can remain temporarily active also 
when the HIF— o has been eliminated, and so drive a consequent oxygen influx. There is hence 
a certain degree of inertia which increments the initial oxygen amount. The effect becomes 
more pronounced as ga gets larger, more HIF-a populating the cell during the transient 
dynamics, before degradation has occurred. However, provided is small enough, the cell 
reaches a steady state configuration, with no residual amount of HIF-a molecules. For this 
reason, we refer to this configuration as to the normoxia regime. 



Vo,i = V,a = = Wi = E,,i = 0. 



(20) 



Wa,i = no 
Oi + Eos + Hp,i = 1-710 



(21) 
(22) 
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As opposed to normoxia condition, a trivial solution of the mean field prob- 
lem is found when: 



Wa,3 = Vo,3 = V,,3 = = £^.3 = ©3 = (23) 

with in addition W3 = uq and £0,3 + "^^,3 = 1 — uq. In this case, the cell is 
dead: the oxygen amount is zero and the hydroxylases cycle is not active, so 
denouncing the lethargic state of the cell device. 

An intermediate asymptotic stationary state is also possible, when the cell 
balances the oxygen consumption by properly adjusting the VEGF production 
cycle, at equilibrium. Upon manipulation of the fixed point equations associated 
to the system (|lip - (|18p . the equilibrium solution for £0,2 is found to verify the 
following quadratic equation: 

T^«-(|-0'^--<|-0^» (.4, 

The previous equation admits a solution only if ,9q, > .g^, otherwise is reduced 
to a sum of three positive terms which cannot mutually cancel for positive co- 
efficients and state variables. When — g^^ then £0.2 ~ Oi which corresponds 
to having the largest allowed oxygen concentration O2 = 1 — tt-q. This lat- 
ter can alternatively be seen as a limiting normoxia condition. All externally 
available cases arc simultaneously occupied by oxygen molecules, which can still 
effectively contrast the HIF-a generation, while providing the necessary supply 
for the cell functioning (Wa,! 7^ 0). Beyond this condition, a quota of HIF-a 
survives the degradation and sustains the VEGF production process. VEGF 
molecules transported outside the cell support the afflux of oxygen and so con- 
tribute to stabilize the quota of HIF-a. In practical terms, the condition gey = gp 
(when the HIF— a and HIE—/? production rates equalize) signals a transition be- 
tween two distinct asymptotic regimes, a fact that could be further appreciated 
by quantifying the values assumed by the remaining involved concentrations. 
We identify the intermediate fixed point with the hypoxia conditions: The cell 
is in fact alive thanks to a dedicated feedback that integrates the original oxygen 
amount via the VEGF production. Once g^ > gp the solution of Eq. (|24)) reads: 



d 

£q.2 = — 



'^fegfi 



(25) 



The equilibrium values of the other variables cannot be univocally deter- 
mined. However, given the value of £0,2 as estimated above, the following ratios 
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are fixed; 



02^2 
Vo,2 
>Va,2^c.2 
Vo,2 

Vo,2 
i£'i,2 

Vo,2 

Notice that an implicit relation between Wa,2 and Vo,2 is also derived by 
making use of Eq. (fT9|) . More interestingly, for the case ga > gp, £0,2 can be 
shown to increase with ga, while keeping all the other parameters constant. The 
maximum value that can be eventually attained is Eo^max = 1 — "fio, a limiting 
condition that, together with ga > gp, identifies the portion of parameters 
plane for which the hypoxia fixed point do exist. More specifically, the largest 
the values of /, i.e. the rate of oxygen production via the VEGF cycle, the wider 
the window in ga that delimits the hypoxia regime. Working in the {ga, /) plane 
we can analytically determine the boundaries of the domain of existence of the 
non trivial fixed point. Its equation can be deduced by imposing £q 2 = £o,max 
into Eq. (|24p which immediately yields to 

.max + (T ed£o .max I 

^ 0«pf2 9a -^2 

This is a line which crosses the axis / = when ga = gp, see Fig. [2] where 
three domains, respectively I, II, III, are outlined. Notice that the normoxia 
and death solutions can in principle extend over all the parameters plane, while 
the hypoxia is solely bound to region II. The hypoxia region shrinks as / is 
reduced and eventually fades off when the VEGF's ability to carry oxygen is 
hindered. We recall in fact that within our simplified formulation, the oxygen 
income is ultimately controlled by the parameter /, and consequently linked to 
the available VEGFq quota. It is hence plausible that that region here identified 
to correspond to the hypoxia regime vanishes when / is set to zero. Numerical 
simulations indicate that the basin of attraction of the normoxia state prevails 
in region I, while the death configuration dynamically takes over in region III. 
In region II, it is the hypoxia condition to represent the stable attractor of the 
dynamics. For this reason, region I, II and III are respectively termed normoxia, 
hypoxia and death, and the caption of Fig. [2] reficcts such a choice. To gain 
further insight into this empirical observation, we turn in the next section to 
studying the stability properties of the identified fixed point. The analysis is 
carried out via combined numerical and analytical means, as the complexity of 
the model prevents us from performing a complete analytical treatment. 
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Figure 2: The existence of the normoxia (I), hypoxia (II) and death (III) regions is displayed 
in the (ga , /) plane, as predicted by the mean-field calculations for N = 1000, no = 20, 
e = 0.9, d = 0.4 and g/j = 0.5. The transition from I to II is found to occur when 

ga — 98 ^ i.e. 

when the HIF— o production rate prevails over the HIF— /3 one. Then the exceeding HIF— o 
is stabilized to an asymptotic fixed concentration via the additional VEGF pathway. It can 
be further appreciated that region II gets compressed when reducing the control parameter 



4-2. Stability analysis 

In this section we set down to study the stabihty properties of the three 
famihcs of cquiUbrium points, as identified above. To this end we calculate 
the Jacobian matrix associated to the system (fTT|) -([T8 l) . explicitly given in the 
Appendix and evaluate the eigenvalues relative to the three fixed points, namely 
the normoxia, the hypoxia and the death. 

Consider first the solution that corresponds to the cellular death, previously 
labeled with the index 3 in the concentration. We recall that Wa.s = V0.3 = 
Vi_3 = Ha.s = £i,3 ~ O3 — 0, Ws — uq, while the concentration £0,3 and 
T-LjS^s can in principle take any value, constrained by the relation £"0,3 + 7^/3,3 = 
1 — no, which ultimately specifies the family of possible solutions. The associated 
characteristic polynomial reads: 



where A stands for the eigenvalue. It is evident that Ai = A2 = A3 = 0, 
A4 = —auQ and A5 ~ —d are possible solution of Eq. (j26p . The remaining 
three eigenvalues cannot be analytically determined. For > gp (a condition 
that embraces region II and III of the parameters plane depicted in Fig. [2]), 
analytical progress is possible: just one stable solution exists among those of 
type 3 that are in principle admissible. To substantiate this observation, let us 



/■ 
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denote with J|"'' the non trivial sub-matrix of the full Jacobian (P^ : 

(-cHp.s +.g„ 

+cMp^3 -e£o,3-d 
+cHp,3 +e£o,3 + d -ga - 9/3 

It is easy to see that det(J|"'') = c'Hi3{9a — g^)(c?-|-efo,3), and so det(J|"'') > 
VH^^s > and 9a > gp- Positive determinant implies unstable solutions, as can 
be immediately appreciated with simple algebraic considerations. This means 
that the fixed point corresponding to the death state is always unstable when 
9a > 9p, and H/j^s ^ 0, this latter quantity being positively defined. At variance, 
stable solution can manifest in the region specified by condition > 9/3, when 
^^(3,3 = 0- This latter request equivalently implies £0.3 = 1 ~ i^o^ and removes 
the degeneracy (at least in region I and II) by uniquely characterizing the fixed 
point of type III selected by the dynamics. The only stable fixed point of type 
"death" in region II and III, is now completely specified and the three lacking 
eigenvalues can be readily calculated, resulting in Ag = 0, A7 = ~e{l — no) — d 
and As = —9a ~ 9/3- Solutions of type 3 are clearly stable also in region I, 
suggesting that the cell can always meet the conditions that drive its death, an 
asymptotic fate which is ultimately determined by the selected initial datum 
(corresponding to the environmental conditions) . 

However, direct numerical inspection, suggests that within region II, the 
system tends to converge to the hypoxia solution, this latter being hence dy- 
namically favored over the stable death regime. Tracing the respective basins 
of attraction can be tackled via a dedicated campaign of numerical simulations. 
Biologically, it can be hence speculated that for sufficiently large values of the 
parameter the cell opposes the death derive, by compressing its respective 
basin of attraction upon activation of the hypoxia cycle. 

In region I, when ga < gp, an infinite family of solutions of type 3 can 
in principle manifest, depending on the specific initial datum. The associated 
eigenvalues can be calculated numerically by tuning both 7^/3,3 and £0,3 within 
the allowed range of variation, pointing to the stability of the solutions. Again, 
the cell can always die in response to particularly stressing conditions. To com- 
plete the description of region I, one needs to address the stability properties of 
the normoxia solution. In this case, repeating the above strategy, the charac- 
teristic polynomial turns out to be 

/ A + 5Wa,i+cH^.i \ 

X^{\+aOi){X+d)dct \ -cHf3,i \ + e£o.i + d =0 

\ -bWa.i - cHp.i -eSo.i - d X + ga+gp / 

which immediately gives Ai = A2 = A3 = 0, A4 = —aOi and A5 = —d. In this 
case the determinant of the residual 3x3 sub-matrix does not return any sen- 
sible information on the sign of the real part of the three unknown eigenvalues. 
To overcome this limitation we proceed as follows: we first integrate numerically 
the governing system of differential equations and then select the attained equi- 
librium fixed point, as an entry to calculate the full set of eigenvalues. Among 
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the missing eigenvalues, two (which can be complex or real depending on ga) 
always display a negative real part, while the third one is real and can change 
the sign as ga is increased. This situation is illustrated in the inset of Fig. 21 
the real eigenvalue is negative over a range of parameters that practically corre- 
sponds to region I, and it approaches zero for g^ ^ gp- This result holds true for 
different parameters choice (i.e. the value of /) and selected initial conditions. 
It can hence be argued that the normoxia condition is stable within region I and 
that it gets unstable when the boundary of region II, the domain of existence 
of the hypoxia solution, is (continuously) approached. Indeed, within region I, 
the competition between the two attractors of type 1 (normoxia) and 3 (death) 
is dynamically resolved to favor the former, as confirmed by direct numerical 
inspection. 

In conclusion, the proposed model can alternatively yield to three distinct 
dynamical regimes, depending on the specific choice of the chemical parameters. 
The normoxia condition is stable in region I and loses its stability when the edges 
of the domain II are approached. In this latter region it is the hypoxia solution 
to take over, by prevailing over the death regime which instead dominates in 
region III. Interestingly, death is always an attractive state of the dynamics, 
irrespectively of the selected parameters, a tendency that the cell efficiently 
opposes in region I and II. Formally, it is hence not accurate to talk about a 
transition among the three states, the landscape of possible stable attractors 
being more complex as outlined above. However, to keep the message simple, 
we shall hereafter reflect the practical segregation into three distinct zones, as 
virtual transitions. Notice that a punctual mutation of the reference parameters 
could substantially alter the cell functioning and determine its death. In the next 
subsection we turn to numerical simulations to validate the proposed scenario. 
The simulations here performed refer both to the idealized mean-field dynamics 
and to the stochastic scenario which respects the inherent discreteness of the 
chemical medium. 

^.3. Numerical simulations vs. mean field theory 

In this section we report about a campaign of simulations aimed at vali- 
dating the theoretical picture depicted above. Mean-field ordinary differential 
equations (fTT|) - (|T8)) can be straightforwardly integrated via a standard Runge- 
Kutta scheme. In Fig. [3] the time evolution of VEGFi and active hydroxylases 
is reported (solid lines) for different values of the parameters. Panels (a) and 
(b) refer to a choice of the parameters that yield to the normoxia condition: 
Once the transient has died out, the system achieves its asymptotic state where 
VEGF molecules are absent. At variance, Wq attains its equilibrium concen- 
tration, different from zero. The same applies to oxygen (not displayed here). 
When increasing the parameter g^^ which controls the HIF— a production rate, 
a transition occurs towards what we termed the hypoxia regime, see panels (c) 
and (d). Here the HIF— a quota is sensibly different from zero, implying the ex- 
istence of a positive feedback reaction on the oxygen production, via the VEGF 
molecules. By further augmenting one gradually approaches the boundary 
of the third allowed region, as previously outlined: When ga is large enough, all 
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chemical concentrations drop to zero, including the Wa amount, which is here 
assumed to signal the cell activity. This scenario is clearly displayed in panels 
(e) and (f) of Fig. O 

Beside integrating the system of differential equations pT|) - ([T5|) . one can 
turn to direct stochastic simulations, which enables one to treat the interacting 
molecules as diffusing massive entities, so respecting the finiteness of the chem- 
ical environment. To this end, we resort to the Gillespie algorithm: as already 
recalled, this is a sophisticated Monte Carlo scheme which enables to inspect 
the actual dynamics of an ensemble made of discrete actors, whose relative in- 
teractions are specified by an assigned set of chemical rules. The outcome of 
the stochastic integration is superposed in Fig. [31 to the curves determined 
within the mean-field approximation. The intrinsic stochastic noise material- 
izes as a local disturbance of the reference mean-field profile, suggesting that 
the discreteness of the environment contributes modestly to the system dynam- 
ics. As opposed to this vision, we will prove in the following that finite size 
correction can be important and can eventually result in drastic modification of 
the dynamical behavior, as predicted within the mean-field scenario. 

Before turning to analyze these aspects, we provide in Fig. |4]a global picture 
of the distinct dynamical regimes that can be faced when tuning the {f,ga) 
parameters: Here the asymptotic concentration of VEGF^, as calculated upon 
integration of Eqs. PT|) - is plotted with a colorcode. The existence of 
three different regions, is clearly confirmed, the transition from normoxia (I) 
to hypoxia (II) being sharper and hence more evidently appreciated. As a 
further complement, we enclose a section of the scanned parameters plane, by 
reporting the recorded concentration of VEGFi versus ga, at fixed /. The 
double transition is shown, bearing intriguing analogy with the phenomenon of 
re-entrant phase transitions. 

5. Dynamics in a multicellular environment: Synchronization and 
cooperative effects 

To shed light onto the cooperative effects that might eventually arise in a 
colony of interacting cell units, we here consider A^c = 3 cells obeying to the 
chemical equations hypothesized above. Now the oxygen fills a shared reservoir, 
from which any individual cell can draw. The cells are hence mutually, though 
indirectly, communicating via the oxygen amount. A schematic layout of the 
interaction network is depicted in Fig. [5l In this case we need to introduce an 
extra label to specify the cell to which the molecules belong to. In practical 
terms, one is forced to deal with the molecules VEGF:^ . Fif , HIF-a^ and HIE-/?-' 
with j = 1, . . . , Nc- The governing chemical equations follow as a straightfor- 
ward generalization of Eqs. H])-®. The mean-field system calculated from the 
underlying master equation is now composed by 4 x A^^ + 3 equations (not ex- 
plicitly given here) , the coupling being indirectly provided by the species hosted 
in the outer environment which constitutes a shared reservoir. 

Working within this generalized setting, interesting cooperative effects arise, 
as supported by the discrete component of the dynamics. In the finite TV regime. 
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correlations exist between the evolving species and translate into important con- 
tributions that can affect the system dynamics. The correlations are instead ne- 
glected within the mean-field approximation, formally valid when the container 
is filled with an infinite amount of molecular constituents. To demonstrate the 
importance of finite size corrections, we assign to two cells a set of parameters 
(including the ga value) which is found to yield to the death condition, in the 
mean field regime and according to the preceding classification. The third cell is 
instead imagined to bear a ga parameter in the hypoxia domain (while keeping 
the other parameters unchanged), as identified above. What is going to happen 
to the considered system of three cells? 

In Fig. ini the evolution of hydroxylase is reported as a function of time. 
The solid lines (red online) stand for the mean-field dynamics, obtained by 
numerical integration of the generalized set of ordinary differential equations 
which accounts for the three inspected cells. Two cells are destined to die (no 
residual Wa), while the third survives (Wa are eventually present). Interestingly 
the stochastic dynamics is peculiarly different, as displayed in Fig. [S] (wiggling 
curves, cyan online). As resulting from a cooperative mechanism, ultimately 
driven by the finiteness of the simulated medium, and so not captured within 
the corresponding continuous description, the two cells destined to death are 
rescued by the third cell. This latter can sustain the angiogenesis process and 
so integrate the oxygen quota which is necessary for the correct functioning of 
its nearest neighbors. The rescued cells do partially contribute by stimulating 
the production of VEGF. Synchronized oscillations are observed in the recorded 
dynamics of the active hydroxylases, as displayed in the zoomed inset of Fig. [6l 

In conclusion, and based on the above findings, it is tempting to speculate 
that the observed degree of cooperation turns out particularly crucial for un- 
derstanding the adaptation of tumoral cells in a changing environment. Under 
stressing condition, it is well known that tumoral cells can undergo mutation, 
changing their specific abilities and fully exploiting their inherent fiexibility. 
Assume that under stressing condition the oxygen income is not enough to 
guarantee the cell functioning. Then the tumoral cells would undergo a criti- 
cal phase, which could induce punctual mutations. In particular, imagine that 
following such random adjustments, a cell enters the regime of hypoxia (in our 
model, by properly tuning the corresponding ga parameter). This corresponds 
to reducing the HIF— a production rate, so consequently freeing a quota of oxy- 
gen, otherwise employed in the HIF— a degradation pathway. Therefore not 
only the tumoral cell interested by the mutation can survive, but also the cells 
populating the local neighborhood (which are plausibly also belonging to the 
tumoral family although not necessarily affected by the same mutation), fol- 
lows a cooperative interaction. This effect is here reproduced within a sound 
dynamical picture where finite size effects prove to be crucial. 

6. Conclusions 

In this paper we have discussed a simplified model to quantitatively address 
the dynamical mechanisms which ultimately underlie the angiogenesis process. 
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The dynamical interplay between a limited number of molecular species is re- 
solved and investigated both via numerical and analytical techniques. The level 
of complexity of the proposed description results from a trade-off balance be- 
tween the request of incorporating a fair amount of biological information, and 
the need of making the problem analytically tractable. Within the proposed for- 
mulation, the cell can follow different routes, depending on the specific values 
of the involved chemical parameters. Distinct macroscopic regimes are here dis- 
cussed and associated to normoxia, hypoxia and death conditions respectively. 
The rate of HIF— a production, as well as the parameter that controls the ability 
of VEGF to call for new oxygen, are shown to play a crucial role in determining 
the transition between the admissible dynamical regimes. These latter are an- 
alytically explained via a mean-field based calculation, which is fully confirmed 
by direct numerical simulations. The role of stochastic noise is also elucidated 
and shown to drive a degree of synchronization between the cells belonging to 
a given colony. This is a spontaneously emerging feature, stemming from the 
intimate discrete nature of the molecular environment, which can potentially 
relate to the attested adaptivity of localized families of tumoral cells. In fact, it 
is in general believed that cancer cells can react to the external stressing stimuli 
(as for the lack of oxygen) by favoring a specific mutation, which renders them 
more adaptable to the hosting environment. While this is a customarily evoked 
scenario, the simultaneous occurrence of the same mutation in all the cancer 
cell population appears a highly improbable event. Cooperative survival mech- 
anisms, as those here identified, could in turn contribute to explain the observed 
robustness, via a self-consistent dynamical argument. As an additional point, 
the results of this paper can be of interest in studying the growth of a tumor 
cord characterized by an aerobic-anaerobic switch triggered by hypoxia [H. 

7. Acknowledgments 

The work was supported by grants from the Associazione Genitori contro le 
Leucemie e Tumori Infantili Noi per Voi, Associazione Italiana per la Ricerca sul 
Cancro, Istituto Toscano Tumori to AA, Ente Cassa di Risparmio di Firenze to 
the Dipartimcnto di Patologia c Oncologia Spcrimcntali, Univcrsita dcgli Studi 
di Firenze. 



Appendix A. Jacobian matrix 



J 



The Jacobian matrix associated to the system ([TT])- 


(fT8|) reads 








/ -aW 
















f£o 







~aW 


-aO 




















+aW 


+aO 


-bHa 


-bWa 
























-bWa - C-Hp 


-cHa 








ga 













-cHp 


-cHa 








.9/3 













+cHf, 




-efo — d 



























-d 







I 





bUa 


bWa + cUp 


cMa 


e£o + d 





-ga - gp 



16 



References 



[1] Astanin, S., Prcziosi, L., 2009. Mathematical modelling of the warburg 
effect in tumour cords. J. Theor. Biol. 258, 578-590. 

[2] Dauxois, T., Di Patti, F., McKane, A., FanelH, D., 2009. Enhanced stochas- 
tic oscillations in autocatalytic reactions. Phys. Rev. E 79, 036112. 

[3] Di Patti, F., Fanelli, D., 2009. Can a microscopic stochastic model explain 
the emergence of pain cycles inpatients? J. Stat. Mech. 1, P01004. 

[4] Fang, J. S., Gillies, R. D., Gatenby, R. A., 2008. Adaptation to hypoxia 
and acidosis in carcinogenesis and tumor progression. Seminars in cancer 
biology 18, 330-337. 

[5] Gardiner, C. W., 1985. Handbook of Stochastic Methods, 2nd Edition. 
Springer. 

[6] Gillespie, D. T., 1976. A general method for numerically simulating the 
stochastic time evolution of coupled chemical reactions. J. Comp. Phys. 22, 
403-434. 

[7] Koong, A. C., Mehta, V. K., Le, Q. T., Fisher, G. A., Terris, D. J., Brown, 
J. M., Bastidas, A. J., Vierra, M., 2000. Pancreatic tumors show high levels 
of hypoxia. Int. J. Radiation Oncology Biol. Phys. 48, 919-922. 

[8] McKane, A. J., Nagy, J. D., Newman, T. J., Stefanini, M. O., 2007. Ampli- 
fied biochemical oscillations in cellular systems. J. Stat. Phys. 128, 165-191. 

[9] McKane, A. J., Newman, T. J., 2005. Predator-prey cycles from resonant 
amplification of demographic stochasticity. Phys. Rev. Lett. 94, 218102. 

[10] Semenza, G. L., 2003. Targeting HIF-1 for cancer therapy. Nature reviews 
3, 721-732. 

[11] Tatum, J. L., Kelloff, G. J., Gillies, R. J., 2006. Hypoxia: importance 
in tumor biology, noninvasive measurement by imaging, and value of its 
measurement in the management of cancer therapy. Int. J. Radiat. Biol. 
82, 699-757. 

[12] van Kampen, N. G., 1992. Stochastic preocesses in Physics and Chemistry. 
North Holland, Amsterdam. 

[13] Vaupel, P., 2004. The role of hypoxia-induced factors in tumor progression. 
The Oncologist 9, 10-17. 



17 




Figure 3: Simulated dynamics of system Ijllll - l|18|l for different choices of the chemical pa- 
rameters. Panels (a) and (b) report the time evolution of the concentration of VEGFi and 
Wa respectively, for a=0.01, b=0.01, c=0.05, e=0.01, /=0.0005, <i=0.0001, g„=0.000002, 
=0.0008. This is the normoxia condition, when the cell functions, beyond the transient, 
without activating the hypoxia cycle. In panel (c) and (d) the homologous quantities (resp. 
VEGFi and Wa) are depicted for Qa = 0.0017, while keeping the other parameters unchanged. 
A residual quota of VEGF; molecules is present, which sustains the hypoxia pathway as de- 
scribed in the main text. Finally, when further increasing (go, = 0.017 in panels (e) and (f), 
respectively referring to the VEGFi and Wa concentrations), the cell dies, no activity being 
asymptotically detected. Solid lines refer to the integration of the mean-field dynamics l|ll|l - 
mSI I. while the associated irregular profiles stand for the corresponding stochastic simulations. 
For all the plots N = 10000, the initial number of O2, W, W* and HIF-a is 1000, that of 
HIF-/3, VEGFi, VEGFo and Ei is equal to 100 and, hence, the initial number of Eq is 5600. 
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Figure 4: Left panel: The asymptotic concentration of VEGF^ is plotted for different values of 
(Sqi/)- The data follows from a direct integration of the mean-field Eqs. |(TTJ - l|18p and are 
represented using the colorcode displayed in the legend. Right panel: The concentration of 
VEGFi is reported as a function of ga, for / = 0.005 (the other parameters are set as in Fig. |3] 
A double transition is clearly displayed, the VEGFi acting as a stabilizer of the cell dynamics 
for a compact range of ga values. Inset a non trivial (see text) real eigenvalue associated 
to the first family of equilibrium points (normoxia), as a function of ga- The eigenvalue is 
seen to approach zero, i.e. the threshold of stability point, when ga = gp- The eigenvalue is 
calculated via an ad hoc semi-analytical procedure which is described in the main body of the 
paper. 




Figure 5: Cartoon of the interaction scheme of three cells. Oxygen molecules (circles laying 
outside the cells' walls) are now shared, and provides an indirect couplings among the three 
cells. Each cell is then supporting a whole cycle of internal reactions, as already exemplified 
in Fig. [T] For the symbols legend, see caption of Fig. [l] 



19 



0.04 
0.035 

0.03 
0.025 

0.02 
0.015 

0.01 
0.005 


-0.005. 



1 1 


1 1 

xlO"^ 


6 
4 






2 















6.5 7 7.5 






1 1 1 1 



4 6 
Time 



X 10 



10 

6 



Figure 6; Wa evolution as a function of time, for a system of three cells, sharing the same 
oxygen reservoir. The parameters are set so to have two cells in the mean-field region deputed 
to death, while the third cell is made to function in the hypoxia domain, via a appropriate 
mutation of the parameter Qa- All other parameters are identical for the three cells. Solid 
lines (red online) represent the moan-field dynamics, while the irregular curves (cyan online) 
stand for the stochastic simulations. Inset: a zoom of the dynamical evolution is represented. 
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